Self-limited self-assembly of chiral filaments 
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The assembly of filamentous bundles with controlled diameters is common in biological systems and desir- 
able for the development of nanomaterials. We discuss dynamical simulations and free energy calculations on 
patchy spheres with chiral pair interactions that spontaneously assemble into filamentous bundles. The chirality 
frustrates long-range crystal order by introducing twist between interacting subunits. For some ranges of system 
parameters this constraint leads to bundles with a finite diameter as the equilibrium state, and in other cases 
frustration is relieved by the formation of defects. While some self-limited structures can be modeled as twisted 
filaments arranged with local hexagonal symmetry, other structures are surprising in their complexity. 



Filamentous bundles assembled from protein subunits are 
essential structural and regulatory components of cells and tis- 
sues. For example, filamentous actin, microtubules, and inter- 
mediate filaments assemble and disassemble to create a strong 
but dynamic cytoskeleton, fibrogen subunits assemble into fib- 
rin fibers and networks to form blood clots (e.g. 01 01) and 
sickle hemoglobin assembles into fibers that impair red blood 
cell function (e.g. JU-JU]). In vivo and in vitro studies sug- 
gest that fibers with finite diameters are the stable morphology 
for fibrin 01 0] under a variety of conditions, while 7-double- 
strand bundles of sickle hemoglobin are metastable JH 0, HI] , 
but the forces that limit filament growth remain unclear. Theo- 
retical calculations ifM I9 UT2I1 have suggested that finite bundle 
diameters are the thermodynamically favored state for twisted 
bundles assembled from chiral subunits. These calculations, 
however, assume specific packings of protofilaments without 
defects, which could relieve strain and thereby enable un- 
bounded growth. The objective of this article is to determine, 
without assumptions about assembly pathways or assemblage 
geometries, if chirality can result in stable bundles with finite 
diameters. We construct a model subunit with simple pair- 
wise chiral interactions that drive assembly into filamentous 
bundles, and combine umbrella sampling lfl3ll and forward 
flux sampling Hill to explore the structures that spontaneously 
assemble for varying degrees of chirality. The simulations 
demonstrate that chirality can result in regular self-limited 
bundles for a range of interaction strengths, but that stronger 
interactions enable defects which give rise to branched net- 
works or irregular bundles. 

Drawing conclusions about self-limited growth in a macro- 
scopic system from simulations with a finite number of sub- 
units is challenging-one must distinguish between simula- 
tions in which growth terminates due to physical constraints 
from those in which the system runs out of subunits lfl5ll . To 
overcome this limitation, we simulate the grand canonical en- 
semble (pVT), in which growth cannot terminate because of 
subunit depletion, since the system is coupled to an unlimited 
bath of free subunits. Before discussing the simulations, we 
consider the conditions for self-limited filamentous assembly 
at fixed total subunit concentration (the NVT ensemble). 

We build on the theory for cylindrical micelles O^HUl] to 



write the free energy for a twisted bundle comprised of rif fil- 
aments and n subunits as G(n, rif) = n(g ro d(nf) — G(l, 1)) + 
G cap (rt, rif) with g ro d the subunit chemical potential in the 
cylindrical 'body' of the bundle and G cap the excess chemical 
potential of the subunits in the 'caps' at either end. G cap and 
g m d are independent of length for long bundles, but depend 
on the bundle diameter, or rif, since subunits at different radii 
experience different environments. For a solution with fixed 
total subunit density ct, minimizing the solution free energy 
density G tot = J2n n, P( n > n f)G(n, rif) — TS with the mixing 
entropy S = -k B J2 n ,n r Pn,m ln(p(n, n f )a 3 ) with a the sub- 
unit size gives the law of mass action result for equilibrium 
bundle densities 



p(n, rif) = exp\j3(pn - G(n, rif))] 



(1) 



with j3 = \/k^T and the monomer chemical potential p = 
G(l, 1) + k&T\n(p(l, l)c 3 ). The onset of spontaneous as- 
sembly is identified with the concentration at which half of 
the subunits are assembled into bundles with the other half 
free in solution, which is given by by lfl8ll c csc sa exp[/3(<7 ro d — 
G(l, 1))]. By minimizing G to t with respect to n and rif we ob- 
tain that bundles follow an exponential distribution of mass- 
averaged lengths P(l) ~ Ze~'/ W, with I = n/rif and the mean 
length (I) ~ [ct exp(Gcap)] 1 ' 2 . However, if there is a bun- 
dle diameter (rif ) that minimizes g ro d, the distribution will be 
sharply peaked about rif when (I) ^> 1, and growth will be 
self-limited. Thus, we begin by measuring g ro d as a function 
of rif. 

We note that in Ref. Il5ll filamentous bundle formation 
is described in terms of the theory of linear polymerization 
H9II20I1 and bundling of linear polymers. That analysis would 
be complicated for our model since the subunit binding free 
energy depends crucially on the number of bundled filaments 
due to the twist imposed by chirality. 

Subunit Model. We consider spherical subunits of diame- 
ter a that are endowed with a polar orientation unit vector p 
and have azimuthal symmetry, ( Fig. [I]). Subunits have pair- 
wise interactions that drive north-to-south pole alignment into 
filaments and weak equator-to-equator interactions that drive 
bundling of filaments. The interaction between two subunits i 
and j is given by 
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V(i,j) = V h ( rij ) + Fp(cos0,dii) 

+V P (cos 9, dji ) + V e (fij , pi , pj ) 



(2) 



2 



where fij = fj — f is the interparticle displacement, nj = 
\fij\, cos 6 = pi ■ pj gives the angle between the two polar 
directions, and dij = \(fj — apj/2) — (f + api/2)\ is the 
distance between poles. Excluded volume is imposed by a 
hard-sphere interaction 
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The pole-pole interaction is given by 

F p (cosM) = -f-pH(d p - d)f{e,e m , x ) 



(3) 



(4) 



with H(x) the Heaviside step function, e p the pole-pole inter- 
action strength, dp the distance tolerance, and max the angle 
tolerance. Parallel alignment is driven by 



f{y,yn 



fexp(-2/ 2 /ymax) V < 2/max 

\ y > 2/ max 



(5) 



The equatorial interaction is given by 



V e (fij,pi,pj) = -e e H(d e - r.y)iJ(/3 max - ,%) 

H(j3 max - @ji)f((p - <y5skew, ¥>max) (6) 

where e e is the equatorial interaction strength, d e is the equa- 
torial distance tolerance, and the width of the equatorial inter- 
action band is set by /3 max with cos/% = |fij • pi\/rij. The 
final factor in Eq. [6] measures the degree of twist, with tp as 
the dihedral angle between the plane (fij, pi) and the plane 
(fj,Pj) (Fig. [T}, which can be calculated from the following 
relations, with qij = fj x pi 



sin tp 



(qij x qji) ■ fi 
Wij\\qji\\rij\ 



cos tp 



Qij ■ Qji 
\qij\Wji\ 



(7) 



The degree of chirality is dictated by the preferred skew an- 
gle t^skew, with <y9 max as the tolerance for deviations from the 
preferred skew. 

Simulations. We explore assembly with Monte Carlo, with 
random translations and rigid body rotations of subunits ac- 
cepted according to the Metropolis criterion lfl3ll . We focus 
on parameters for which nucleation is a rare event and bun- 
dles grow by addition of monomers, and thus we use only 
single particle moves. To ensure that self-limited growth is 
not a result of subunit depletion, we sample the grand canon- 
ical ensemble by coupling the system to a subunit bath at 
chemical potential ^ with subunit insertion/deletion moves 
lfl3tl . Since there are large nucleation barriers, we employ 
umbrella sampling simulations! 13] to calculate the free en- 
ergy. We use the total bundle size n as the reaction coor- 
dinate, and measure the probability p(n, ni) that a particu- 
lar subunit is in a bundle of size n, using a series of win- 
dows in which hard walls constrain the size of that cluster to 
a range of n. The free energy is then obtained from Eq. Q] 
with G(n, n,f) = —ksT \n[p(n, rif)(T 3 } + ^ with p(n,rif) — 
p(n,n f )/nMM- 

The free energy for <y9 s k e w = 0.38 and snapshots of repre- 
sentative bundle configurations are shown in Fig. [2] We see a 




FIG. 1: Subunit model, (left) The pole-pole interaction, dij is the 
distance between the corresponding poles, and 9 is the angle be- 
tween f>i and pj. (right) The equatorial interaction, tp is the di- 
hedral angle between the plane (fj,Pi) and the plane (fij,pj). For 
all simulations reported in this work, the polar interaction strength 
is e p = 14/cbTo, the lateral interaction strength is e e = 2.95/cbTo, 
and the distance and angle tolerance parameters are d p = O.lcr, 



d c = O.lcr, ^ max = 1, 



= #max = 0.25 with angles in ra- 



dians. The temperature is T = To, except for the simulation of 
V?skew = 0, for which T — 1.05To. The GC bath chemical potential 
is (jtb = feBTlnO.01. For this model G(l, 1) = fc B Tln87r 2 due to 
rotational entropy. 



rapid rise in free energy at small n during which short struc- 
tures with rif = 3, 4 , and 5 filaments appear successively, fol- 
lowed by the critical nucleus with rif = 7 and n w 25 subunits 
(Fig. |2j\. The free energy is unfavorable below this size be- 
cause the majority of subunits have unsatisfied lateral and/or 
polar contacts (see Cap Free Energy below). After reaching 
the critical nucleus, the bundle grows lengthwise in both direc- 
tions, while maintaining the same structure, and the free en- 
ergy decreases linearly. As shown below, the Uf = 7 structure 
corresponds to the optimal bundle structure for <£> s k e w = 0.38 
and thus further lateral growth is unfavorable. The slope of 
the free energy in this region corresponds to the chemical po- 
tential S-rod(^f = 7). 

While the exponential distribution of filament lengths is 
derived above for the NVT ensemble, in the /iVT ensemble 
the bundle will continue to grow lengthwise indefinitely (pro- 
vided that g ro d < —fib). To evaluate the free energy of lateral 
growth, we imposed hard spherical boundary conditions with 
a diameter D = 44a, which is large enough that bundle prop- 
erties are independent of D. Upon reaching the boundary, the 
bundle grows by increasing its diameter, which results in the 
increasing free energy at large n in Fig. [2] 

Cap free energy. The cap free energy is calculated using 
Gca P (™, rif) = G(n, ni) - n(g md (n.f) - g x ), with g md obtained 
from the slope of the linear regime in the free energy (or from 
Fig. [3] below). As shown in Fig. [2] G cap rises rapidly un- 
til saturating at the critical nucleus with G cap « 37fcBT for 
Vskew = 0.38. This value is similar to cap energies measured 
for cylindrical micelles ifTill and corresponds to a large aver- 
age bundle length in the canonical ensemble; using the mea- 
sured 5 ro d and G cap (rt, rif) we solved Eq.Q]to obtain a mass- 
averaged bundle size of about 10 8 subunits at the CSC. Al- 
though the magnitude of G cap depends on the strength of the 
polar bonds (e p = 14/cbT), we find that robust bundle forma- 
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FIG. 2: . (left) Free energy G and cap free energy G cilp as functions 
of the number of subunits n with i^skew = 0.38 obtained from um- 
brella sampling. The rise in free energy at large n occurs when the 
bundle reaches the hard system boundary and begins to grow a third 
layer, (right) Structures corresponding to the indicated points on the 
free energy plot. 



tion requires e p ^> e e , implying that large cap free energies 
and hence large average filament lengths are general. 

Similar results were obtained for lower skew angles 
Vskew = 0.32 and <y9 s k ew = 0.25 at low n, but convergence 
in the free energy calculation was questionable because tran- 
sitions between different values of rif were rare at large n. We 
overcame this limitation as follows. 

Self-limited bundle diameters depend on preferred 
skew. In the umbrella sampling simulations the system adopts 
the number of filaments n£ that minimizes G for a given n. To 
determine the dependence of the subunit free energy on 
rif, we performed additional sets of 'constant filament num- 
ber' umbrella sampling (CFNUS) simulations in which rif and 
n are constrained. A small structure with rif filaments is ex- 
tracted from an umbrella sampling simulation, and subjected 
to a simulation in which any move that changes rif is rejected. 
Specifically, moves which cause the number of subunits in any 
filament to differ by more than 3 are rejected and the mean 
of each filament along the bundle axis must remain within 
2(7 of the bundle center. These additional requirements con- 
strain the configurations of the cap and hence affect G cap , but 
do not affect <7 r0 d in long bundles, which is determined from 
Srod(^f) = {dG(n, nf)/dn), h - at large n (when g m ^ becomes 
independent of n). This procedure is repeated for all com- 
monly observed morphologies with a given rif. 

The chemical potential <? ro d calculated from the CFNUS 
simulations is shown as a function of bundle size for three 
preferred skew angles in Fig. [3] In each case, there is an opti- 
mal bundle diameter, or number of filaments rif . Although the 
minima appear shallow, the large average bundle lengths cal- 
culated above ensure that the free energies for different bun- 
dle morphologies differ by many k^T. Note that nj? = 7 
for <£ s kew = 0.38, in agreement with the unrestrained um- 
brella sampling. Furthermore, the slope of the free energy 
in the linear region of Fig. [2]gives g IO d = —4.75 (via Eq. [TJ, 
which matches the chemical potential determined for rif = 7 
in Fig. [3] showing that the two protocols agree. 

The existence of an optimal diameter can be understood as 
follows. Adding layers decreases energy, since subunits in the 
outermost layer have unsatisfied lateral contacts. However, 
the preferred skew (p s k ew causes filaments to tilt with respect 




Bundle structures 



FIG. 3: (left) Subunit free energy as a function of filament struc- 
ture for three skew angles from the constrained umbrella sampling 
(CFNUS). Structures are labeled with the number of filaments in 
each layer starting from the center. Results are shown only for the 
lowest free energy morphology at each value of rif. (right) Snapshots 
from CFNUS simulations, shown in cross-section view, illustrate op- 
timal bundle morphologies for some values of rif. 



to the central filament; tilt increases with layer number for 
preferred skew angles </? s ke W above a critical value. Due to 
tilt, filaments trace a curved path around the bundle, which 
requires unfavorable bending of polar bonds. Furthermore, 
complete formation of lateral bonds between layers, requires 
that the exterior filaments stretch while the interior filaments 
compress. With each additional layer, the degree of exten- 
sion and compression increases. These effects overwhelm the 
energetic benefit of adding an additional layer at the optimal 
bundle diameter nj? . Since the magnitude of the tilt increases 
with ^ s kew; n f decreases with increasing <y9 s k ew (Fig. |3J. In 
agreement with this explanation, the chemical potential for 
<^skew = (corresponding to achiral interactions) does NOT 
show a minimum, and as shown below (Fig. [4b ) bundles with 
<^skew = have unbounded lateral growth. This observation 
confirms that chirality is the reason for self-limited bundle 
sizes in this model. 

Interestingly, the optimal bundle morphology for a given 
number of filaments changes with rif. As shown in Fig.[3]3, the 
central layer of the bundle can vary between 1 and 4 protofil- 
aments, and usually corresponds to the structure that maxi- 
mizes rotational symmetry. While these are the lowest free 
energy structures at each value of rif among those taken from 
unrestrained umbrella simulations, we cannot rule out lower 
free energy morphologies that we did not test. 

Dynamics. Having shown that self-limited bundles are 
the equilibrium state for our model chiral subunits, we now 
demonstrate that they are also the kinetically selected state. 
Since bundle formation is not accessible by straightforward 
dynamical simulations due to the large nucleation barrier 
(Fig. [2), we used forward flux sampling (FFS) ifTill to obtain 
an unbiased ensemble of assembly trajectories. We used the 
bundle size n as the order parameter and performed FFS until 
bundles reached a size larger than the critical nucleus (it is not 
necessary that the order parameter be a good reaction coordi- 
nate, although a bad order parameter can inhibit convergence). 
At this point, FFS was no longer needed, and we continued the 
simulations with straightforward dynamic Monte Carlo; nu- 
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f [sweeps] 

FIG. 4: (a) Bundle size as a function of Monte Carlo (MC) sweeps 
for three trajectories with ^skew = and T = 1.05Tb. Trajectories 
were initiated using forward flux sampling (FFS) as described in the 
text, dynamics are shown after FFS ended. The solid line indicates 
the scaling n(t) ~ t 2 . (b) A snapshot from one trajectory in (a). 




FIG. 5: A branched bundle formed in a dynamical trajectory with 
Vskew = 0.32. 

cleated bundles readily grow in length until they reached the 
imposed boundary (D = 32a), but lateral growth terminates 
at m = 7, in agreement with the equilibrium calculations. 
The dynamic nucleation pathways observed with FFS closely 
follow the minimum free energy pathway observed with um- 
brella sampling (Fig. EJ?), indicating that structures below the 
nucleus size achieve relative equilibrium quickly in compari- 
son to the nucleation time l23l l24ll . 

In contrast, assembly trajectories for <^ s k e w = (Fig. |4]i 
demonstrate unbounded lateral growth, and after reaching a 
size of n « 100 scale as n(t) ~ i 2 , with t the number of 
Monte Carlo sweeps. This scaling is consistent with growth 
dominated by subunit addition to the bundle body, and g ro d 



independent of rif, in agreement with the free energy calcula- 
tions. Although the bundle grows with hexagonal order, we 
note that pentagonal defects become trapped within the as- 
semblage, particularly during rapid growth that occurs under 
stronger interactions. 

Strong interactions lead to branched networks. Bundles 
with 7if > rif tend to form defects that relieve strain. In some 
cases these defects serve as nucleation points for branching, 
which enables further growth until the branch reaches its op- 
timal diameter and a system boundary. As shown in Fig. [5] 
each branch maintains a finite radius. Additional branc hing 
and bundle growth lead to a branched network (see Ref. IU9I1 
for further discussion of branched bundles). We also observe 
branching in the canonical ensemble simulations with strong 
interactions. While this simple model is not intended to rep- 
resent a particular molecule and crowding affects assembly 
at high density JfJ], we note that fibrin clots are composed of 
branched networks of fibrin bundles with uniform bundle radii 
(e.g. H). 

In conclusion, we simulated the assembly of subunits with 
a simple potential that drives the formation of filamentous 
bundles. For moderate interaction strengths, the local pack- 
ing constraints that arise due to chirality cause assembly to 
terminate at a finite bundle diameter, while bundles propa- 
gate easily in length. Stronger interactions, however, lead 
to defects which enable the formation of multiply connected 
branched networks. The optimal morphologies of assem- 
bled bundles with different numbers of filaments have differ- 
ent symmetries. The simulation results indicate that sponta- 
neously assembled structures can deviate significantly from 
regular hexagonal bundles, and thus it is important to evalu- 
ate assembly behavior with dynamical algorithms that do not 
impose particular assembly pathways or morphologies. The 
approach we have adopted to evaluate self-limited growth in a 
finite-sized simulation could be used to understand specific bi- 
ological molecules; for example, a patchy-sphere model could 
be constructed from atomic -resolution structures of sickle 
hemoglobin in order to understand the effects of chirality and 
sphere-packing on hemoglobin filament assembly HHH]. 
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